Vibrationally induced flip motion of a hydroxyl dimer on Cu(llO) 
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Recent low-temperature scanning-tunneling microscopy experiments [T. Kumagai et al, Phys. 
Rev. B 79, 035423 (2009)] observed the vibrationally induced flip motion of a hydroxyl dimer 
(0D)2 on Cu(llO). We propose a model to describe two-level fluctuations and current- volt age 
characteristics of nanoscale systems which undergo vibrationally induced switching. The parameters 
of the model are based on comprehensive density-functional calculations of the system's vibrational 
properties. For the dimer (0D)2 the calculated population of the high and low conductance states, 
the I -V, dl/dV, and d^I/dV^ curves are in good agreement with the experimental results and 
underlines the different roles played by the free and shared OD stretch modes of the dimer. 



Electron transport through single-molecule junctions 
has been receiving enthusiastic interest for a development 
of novel molecular devices. Nonlinear I — V character- 
istics associated with the vibrationally mediated config- 
urational change with different conductances have been 
observed in a series of systems such as pyrrolidine on a 
Cu(001)\ H2 on Cu^, CO bridging a Ft contact^ and H2 
in An contacts^. In these systems dl/dV spectra show 
anomalous spikes — in contrast to steps usually observed 
in inelastic electron tunneling spectroscopy (lETS)^ — at 
the bias voltage related to the vibrational mode energies. 

Recently Kumagai et al.^ studied the dynamics of 
a single hydroxyl (OH, OD) molecule and the dimer 
(00)2 on Cu(llO) using a scanning tunneling microscope 
(STM). The STM images observed for the monomer sug- 
gested the possible quantum tunneling of a hydroxyl 
between two equivalent adsorption configurations on 
Cu(llO), as supported by the density functional theory 
(DFT) calculations of the transition path and rate of the 
flipping of OH on Cu(llO)^. This spontaneous flip mo- 
tion of hydrogen atoms in the monomer is quenched for 
the dimer at low temperatures, but can be induced by 
excitation of the OH/OD stretch mode by tunneling elec- 
trons. Time-averaged measurements of the current show 
a non-linear current (/) increase at the bias voltage (V) 
inducing the transition from the high and low conduc- 
tance states. The appearance of the peak in dl/dV and 
the peak and dip in d^I/dV^ from transitions between 
states with distinct conductances have also been reported 
previously^i"— . 

In this paper a combined use of the DFT-based 
SIESTA^, TranSIESTAS and InelasticaiSiii packages 



permits us to gain insight into the elementary processes 
that induce the flip motion of the asymmetric dimer. The 
extensive DFT calculations provide the ground state ge- 
ometry, vibrational modes, electron-vibration couplings, 
emission rate of vibrations from tunneling electrons, vi- 
brational damping due to electron-hole pair excitation, 
and the high and low conductance. These calculated 
properties allow us to model the population of the high 
and low conductance states as a function of the bias volt- 
age and the nonlinear I — V characteristics for (0D)2 on 
a Cu(llO) surface. The experimental results (relative oc- 
cupation of the high and low conductance state, I — V 
curve, and dl/dV) are nicely reproduced, and the differ- 
ent roles played by the free and shared OD stretch modes 
in the vibrationally mediated configurational flip motion 
are clarified. 

The telegraph switching between high- and low- 
conductance states of (0D)2 is shown schematically in 
Fig- HI In low-temperature STM, no spontaneous switch- 
ing is observed at low bias since the barrier for the re- 
action is substantial^. Keeping the STM tip station- 
ary and increasing bias over approximately 200 meV the 
switching between the two degenerate low-energy con- 
figurations is triggered by phonon emission. However, 
the reaction rate remains rather small until the bias ex- 
ceeds the v{OT>) stretch vibrations (w 300-330 meV). We 
denote the two configurations as the high [H) and low 
(L) conductance states where the tip-dimer distance is 
smaller for the H configuration. 

To investigate the flip motion of the dimer (0D)2 we 
performed calculations on a 6 atom thick Cu 4x3 slab 
using periodic boundary conditions. The tip was mod- 
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FIG. 1: (Color online) Schematic view of the high and low 
conductance states of the dimer (0D)2 on Cu(llO). Yellow, 
blue and red circles are Cu, O, and D atoms, respectively. See 
text for details on the relaxed geometry. 



eled as one protruding Cu atom on the reverse side of 
the slab. All calculations were performed with the PBE 
GGA functional, DZP (SZP) basis set for dimer (bulk 
Cu), 3x4 /c-point sampling, and a 300 Ry mesh cutoff. 
The elastic transport properties were calculated using the 
DFT+NEGF method where 6 additional layers of Cu 
were used to connect the central region to semi-infinite 
metallic leads. Fig. [T] shows the schematic side view of 
the dimer in the H and L conductance states. The cal- 
culated low-bias conductance ratio Gh/Gl = 2.0 closely 
matches the experimental ratio. However, the calculated 
absolute conductances are much larger than in the exper- 
iment since numerical considerations demand the calcula- 
tions to be performed at a small tip-dimer distance. The 
relaxed hydrogen-bonded OD-0 distance is 2.90 A while 
the OD bond length and tilt angle to the surface normal 
are 1.00 (0.98) A and 81° (51°), respectively. The two 
hydroxyl groups, bonded near adjacent bridge sites along 
(001) are oriented according to the optimal configuration 
for hydrogen bonding. The D atom pointing towards the 
adjacent OD molecule form the OD-0 bond (with char- 
acteristic shared OD stretch mode) in contrast to the 
D atom pointing away from the adjacent OD molecule 
(with free OD stretch mode). These nonequivalent con- 
figurations of each OD molecule lift the degeneracies of 
the vibrational energies for the two OD molecules. 

The vibrational modes and frequencies of the 
molecular adsorbates and the corresponding electron- 
phonon couplings are calculated from a finite difference 
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These calculations were repeated with the 



tip scanned over the surface as indicated by the black dots 
in Fig. [21 Table I lists the calculated vibrational modes 
for the H and L configurations. Here z^(OD) labels the 
free OD stretch mode, ;/(0D-0) the shared OD stretch 
mode involving the D atom between two oxygen atoms, 
and T:otxy{z) rotation modes in the surface plane (surface 
normal) . Values for the low-energy modes mainly involv- 
ing motion of the oxygen atoms have been omitted from 
the table. 

In addition to the vibrational energies (?iw), the 
electron- hole pair damping rate jeh and vibration genera- 
tion rate 7om were calculated for each tip position within 
the lowest order expansion (LOE) schem e^°i^^ . In the 
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FIG. 2: (Color online) (a) Emission rate constant A^^^ ^(OD) 
of the free i/{OT)) stretch and (b) A^^^ i/(OD-0) "-"^ shared 
!^(0D-0) stretch mode. The blue and red circle are O and D 
atom, respectively. The black doted grids are the tip position 
fixed at 5.6 A above the oxygen atom and the red arrows 
indicate the vibrational motion. The values of the contour 
lines are indicated in the figures [10'"'/(s eV)] and the exact 
values with the tip positioned at the large green dots are listed 
in Table [B 



low-temperature limit the emission rate is give by 

IcmiV) = Ae,n (le^l - kuj) {\eV\ ~ huj) (1) 
where V is the applied bias, 9 the step function, and 
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A, 



the 



- ;i-Tr[MAiMA2] = f 
emission rate constant written in terms of electron- 
phonon coupling matrix M and partia l sp ectral den- 
sity Ai/2 at the Fermi energy, cf. Ref. uQ. The spec- 
tral density can further be rewritten in terms of a sum 

1 /2 

over scattering states \tpa ) at the Fermi-energy incident 
from substrate / tip (1/2) showing the equivalence of the 
LOE scheme with the Fermi-golden rule, assuming a con- 
stant density of states around the Fermi-energy^i'^ Fur- 
thermore, for weak coupling between tip-adsorbate, the 
emission rate constants scale with the tip-adsorbate cou- 
pling (r) squared (A oc r^). Since the current scales in 
the same way, the emission rate at a given voltage is pro- 
portional to the current. The electron-hole pair damping 
rate 'jeh is insensitive to the position of the STM tip 
since the damping is dominated by the metal surface. 
The emission rate map for the two J'(OD) and z^(OD-O) 



TABLE I: Vibrational mode frequencies, 
stants and their electron-hole pair dampin, 
conductance state. 



emission rate con- 
g rates at H and L 



Mode 


Tiuj [meV] 


Aem [lOV(seV)] 


7eh [lOVsl 


i/(OD) L 


327.2 


0.17 


47.3 


v{OT)) H 


326.6 


4.42 


52.3 


u{OT)-0) L 


301.2 


0.49 


248.0 


v{OT>-0) H 


301.8 


0.35 


249.6 


rot,(OD-0) L 


76.8 


1.19 


250.7 


rot,(OD-0) H 


76.5 


0.67 


268.0 


rot^;,(OD-0) L 


77.0 


0.24 


103.7 


rot,y(OD-0) H 


76.9 


0.10 


76.0 


rot^y(OD) L 


51.8 


0.01 


59.4 


rot,y(OD) H 


52.6 


0.71 


60.7 


rot2(OD) L 


49.9 


0.27 


97.0 


rot3(OD) H 


49.8 


0.79 


86.6 



modes are shown in Fig.[5]where the emission rates were 
calculated by scanning the tip over the surface at a con- 
stant height of 5.6 A above the oxygen atoms, i.e., the 
DFT calculations of e-ph coupling, transmission, emis- 
sion rates, and damping were repeated for geometries 
with the tip displaced relative to the OD dimer. We note 
that the emission rate for the free OD stretch is more lo- 
calized around the high conductance site than the shared 
OD-0 stretch mode which is also evident from Tabic T 
The green points in Fig. [5] correspond to the H and L 
configurations listed in Table ID 

To describe vibrationally induced switching in 
nanoscale systems we propose a simple model which ex- 
presses the current in terms of the occupation n^jj^^ and 
conductance crH{L) of the H and L conductance states. 



I = <7HnH{V)V + cTLnL{V)V, 



(2) 



where 71^(^1,^ is determined as a stationary solution of the 
rate equations. 



driL/dt = r^" 
with the condition nn H 



(3) 
(4) 



n-L 



riLiV) 



— 1. This gives 
nH{V) 



■H ' 

(5) 

where F^^^ is the transition rate from H to L. The ex- 
periments show that the reaction rate at a constant volt- 
age is propotional to the current, i.e., the reaction rate is 
proportional to the emission rate of the vibrations. The 
simplest way to model the transition rates is therefore to 
assume a linear dependence on the vibrational occupa- 
tion cx 7cm /7oh or equivalently the vibrational generation 
rates 7enr^, i.e., 

YH^LiL^H)^y~^ = Fo^(^) + C.F.^^y 



{V), (6) 
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FIG. 3: (Color online) (a) Calculated transition rates 
^H^L{L^H)^y^ and (b) occupations H (red) and L (blue) 
configurations of a dimer (0D)2 on Cu(llO). See text for 
the parameters used to fit the experimental occupations from 
Ref.[g (blue/red dots). The dashed curves are calculated with- 
out the contribution of the free OD stretch mode. 



where a constant rate F, 



ff(L) 




is introduced to model the 



collective effect of many low-energy vibration modes and 
where F^^^^ ^ represents the broadened vibrational gener- 
ation rate— of a distinct mode i given by 



VF(cj - rif (^\alh)7cm,»(ey,c^)dcj. (7) 



In the above equation we use a Gaussian distribution 
function VF(cj, Cph) oc /(2c'ph) characterized by a 
standard deviation cTph. The prefactors Ci describe 
the likelihood of the generated vibrational excitation to 
induce the flip motion, i.e., the probability of trans- 
ferring the motion from the vibration to the reaction 
coordinated-^. 

We obtain a good fit to the experimental OD dimer 
data from Ref. considering only three different vibra- 
tional modes, see Fig. [3] The fitting parameters used 
{Ci,m"l^ ,al^) are (0.2 x 10-^ 327.2/326.6 meV, 20 
meV) for i^(OD), (1.0 x 10^^, 301.2/301.8 meV, 8 meV) 
for t/(OD-0), (1.2 X lO-'^, 76.8/76.5 meV, 19 meV) for 
rot^(OD-O), and F^^-^ = 0.7 x lOVl x lO"" s^^. We 
first note that we need to use large values for the broad- 
ening (Tpjj ~ 10 — 20 meV to fit the experimental data. 
We believe this to be caused by the statistical nature of 
the experimental data collection and do not reflect ther- 
mal or phonon DOS broadeningi^. The fitting constants 
Ci and F^^^ are only determined up to a multiplicative 



factor since the occupations are determined by the ratios 
in Eq. We choose the prefactor of the shared i^(OD- 
O) mode as the reference and set it to 10~^ based on an 
order of magnitude estimate from the experimental flip 

frequencj*^. The values for Tq^^ only determine the low- 
bias occupation where :s> Tq because the system is 
experimentally only observed in the L configuration. The 
reason for this preference of L at low bias might be due 
to the larger current and thus larger vibrational genera- 
tion rate of low-energy phonons in the H configuration. 
However, in the intermediate bias range 200 ^ 300 meV 
we see a slight change in the occupation as shown in 
Fig. EJb). Although we cannot unambiguously assign a 
vibrational mode to this change in occupation from the 
available experimental data, we have chosen to model 
this by the rotz(OD-O) mode. In contrast, there is much 
less latitude in the fitting parameters for the z^(OD) and 
i^(OD-O) modes. We note that the change in occupation 
at 300 mcV do not fit with the free z^(OD) vibrational 
energy (w 330 meV) and clearly indicate that the shared 
i'(OD-O) mode is the main culprit in inducing the flip 
motion. This assignment is supported by the fact that 
the occupation quickly approaches 50/50 which implies 
that the emission rate constants Acm of the high- and 
low-conductance states are of similar magnitude. This is 
clearly not the case for the free ;/(0D) but true for the 
i'(OD-O) mode, see Tab. |T] and Fig. H To underline the 
effects of the shared z^(OD-O) and free i'(OD) stretch 
modes Fig. [3] shows the full modeling (solid line) and 
without the free (/(OD) vibration (dashed line). With- 
out the z^(OD) mode, the occupation rapidly approach 
the ratio of the shared i^(OD-O) emission rate constants 
nL/nn > -^mn/'^cm ~ 0.71. Note that the approx- 

imately equal occupation of the H/L states is a coinci- 
dence and that the ratio of the occupations at high bias is 
simply given by the ratio of emission rate constants and 
the fitting constants Ci. In contrast to the shared i/(OD- 
O) mode, which increases the H state occupation, the 
main effect of the free i^(OD) mode is a slight increase of 
the H ^ L transition rate, Fig.[2{a), and thereby slightly 
decreases the occupation of the H state. 

Using the calculated conductances {<jh — 145.62 
nA/V, (Tl = 78.34 nA/V) scaled by the constant factor 
2.56 X 10~^ to account for the small tip-sample distance 
used in the calculation and the bias dependent population 
[nH{V), TiLiy)] shown in Fig. [3Kb), the I ~V character- 
istics were calculated from Eq. As shown in Fig. 0] 
the calculated time averaged I — V curve [Fig. 4(a)] and 
cPl/dV^ [Fig. 4(b)] show remarkable agreement with the 
experimental results. We note that the lineshape of the 
cPl/dV^ signal is clearly different from inelastic electron 
tunneling spectra which normally only shows a peak or 
dip^'*. In addition, the size of the signal from the vibra- 
tion is much larger than what one normally associates 
with lETS spectra. 

In summary, we have shown that the flip motion be- 
tween high and low conductance conflguration of the OD 
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FIG. 4: (Color online) (a) I — V curve (blue line: calculation, 
red dots: experiment) and (b) dl /dV (green line: calculation) 
and d^I/dV^ (black line: calculation, red dots: experiment). 
The experimental data originate from Ref. @. 



dimer on Cu(llO) is mainly induced by the excitation 
of the hydrogen-bonded shared OD stretch mode. Be- 
cause of the unique asymmetric inclined orientation of 
each hydroxyl, the shared and free OD stretch modes 
have different vibrational frequencies and consequently 
affect the flip motion at different applied bias voltages. 
The relative occupations {uh and ul) of the high and 
low conductance state as a function of the bias voltage 
were nicely reproduced by solving a simple rate equation 
for uh and tt-l in terms of the transition rate between 
the H and L configurations. The calculated high and 
low conductance (cr/f , ctl) and their occupations [nH{V), 
nL{V)] enabled us to obtain the nonlinear I — V curve, 
dl/dV and d^I/dV^ in excellent agreement with the ex- 
perimental results. The presented theoretical analysis 
based on extensive DFT calculations (stable configura- 
tions, vibrational modes including their generation rates 
by tunneling electrons and damping rates) is not limited 
to the specific case of hydroxyl dimers on Cu(llO), but 
can also be applied to other systems which exhibit non- 
linear I — V characteristics arising from the vibrationally 
mediated switching between high and low conductance 
states. 
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